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ABSTRACT 

A weakly ionized plasma in a strong and nonuniform magnetic field 
exhibits an instability analogous to the flute instability in a fully ionized 
plasma. The instability sets in at a critical magnetic field. To study 
the final state of the plasma after the onset of the instability, we numerically 
integrate the plasma equations assuming a certain initial spectrum of small 
disturbances . In the regime we studied, numerical results indicate a final 
steadily oscillating state consisting of a single finite amplitude mode 
together with a time independent modification of the original equilibrium. 

Our numerical results agree with the analytic results obtained by Simon in the 
slightly super-critical regime. As the magnetic field is increased further, 
the wavelength of the final oscillation becomes nonunique. There exists 
a sub interval in the unstable wave band. Final stable oscillation with a 
wavelength in this sub interval can be established if the initial disturbance 
has a sufficiently strong component at the particular wavelength. 
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1. Introduction 

The information one obtains from a linear stability analysis is the 
critical value of some external parameter for the onset of the instability and 
the wavelengths of the modes which begin to grow. After a period of exponential 
growth of the unstable modes to finite amplitudes, the linearized equation 
with the assumption of infinitesimal perturbation is no longer valid. 

To learn the ultimate fate of the excited modes and the final state of 
the plasma, a nonlinear treatment is necessary. In the case of instabilities 
in a weakly ionized plasma obeying simple moment equations, Simon 1 has developed 
a general theory describing the final oscillating state of the plasma. 

The initially stable plasma is assumed to become unstable upon small fractional 
increase A of some external parameter. By expanding with respect to the 
small parameter A, the final amplitudes of the steadily oscillating states, the 
frequency shifts and the time independent change in the equilibrium distribution 
can be determined in terms of A* 

Another possible approach to the nonlinear stability problem is through 
numerical integration of the plasma equations using a computer. One assumes 
a certain initial spectrum of small disturbances and studies its development 
in time to learn the behavior of the plasma upon onset of instability. 

The work of Sato and Tsuda 2 is such a study of the cross- field instability. 

In this paper, using a similar computational approach, we study the 
nonlinear evolution of an instability which causes the convective flow of a 
weakly ionized plasma in a strong and nonuniform magnetic field. Part of 
this problem has been discussed by Kadomtsev 3 in the quasi-linear approximation. 
Nonlinear analytic results describing the final steadily oscillating states 
in the slightly super-critical regime, A « 1, have also been obtained by 
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Simon 4 using the general theory of Ref. 1. However, the computational approach 
has the advantage of not involving the assumption of small super- criticality. 

Thus we can study the behavior of the plasma as the external parameter is 
increased further from its critical value for the onset of the instability. 

This is the purpose of the present investigation in addition to comparing 
the numerical results in the slightly super-critical regime with the analytic 
results of Ref. 4 (hereafter referred to as I). 

The instability we studied is analogous to the flute instability in a 
fully ionized plasma. Consider a weakly ionized plasma filling the space 
between two infinitely long concentric cylinders (cf. Fig. l). Let us assume 
that an azimuthal magnetic field ^ decreasing along the radius as l/r, 
exists in the space between the cylinders and that in the equilibrium state 
all quantities vary in the radial direction only. The oppositely directed 
drifts of electrons and ions in the nonuniform field produce no charge separation 
in the equilibrium state. However, a z-dependent density perturbation will 
produce charge separation and associated electric fields. The resultant 
E X H drifts, coupled with an appropriate equilibrium radial density gradient, 
can cause growth of the initial perturbations. Diffusion due to collisions 
with the background neutrals tends to smooth out these perturbations and 
hence the magnetic field has to exceed a certain critical value for the 
instability to occur. 

In Sec. 2, we write down the basic equations of the system and describe 
the equilibrium state as given in I. We then write the equations governing 
the perturbations and give a brief account of the linear stability analysis. 

The perturbation equations are properly scaled in Sec. 3. We then expand 
the perturbations in Fourier series and arrive at an initial value problem 
involving an infinite number of interacting modes. In carrying out the 
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numerical integration on an IBM 360 computer, we approximate the system by 
one with a finite number of modes and consider initial white-noise- like or 
weighted small perturbations. The time development of these perturbations 
at various magnetic field strengths above the critical field is presented 
in Sec. 4 and Sec. 5* We give a final discussion in Sec. 6. In the regime 
we studied, our computed results indicate a final steadily oscillating state 
with the dominance of a single finite amplitude mode together with a time 
independent modification of the equilibrium distribution. The final oscillating 
state approaches that given in I when the magnetic field is slightly above 
its critical value. 

2. Basic Equations and Linear Stability Analysis 

We consider a weakly ionized plasma filling the space between two 
conducting grounded long cylinders with radii R and R + d respectively, where d « R 
(cf. Fig. l). We will assume quasineutrality, i.e. the number densities of 
electrons and ions are equal. Ionization at the two cylinders is achieved 
in such a way that a constant density s is maintained at the inner cylinder 
and a density s - 6s at the outer cylinder, where 6s « s. There is also 
an applied strong steady toroidal magnetic field produced by a current flowing 
along the inner cylinder. If we assume the plasma current to be negligible, 
by Maxwell's equations the magnetic field H is in the 9 direction with ! .' 

intensity decreasing as l/r. 

Since the density of the charged particles is considerably lower than 
the density of neutral particles, we can assume that the neutral component 
is stationary. We shall describe the behavior of electrons and ions by the 
equation of continuity and the equation of motion of each species. They 


have the form 



% + V (nyf)= 0 


T + Vn = + en7<p ± ^ ( v* H ) - — ^ 
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( 2 . 1 ) 


( 2 . 2 ) 


where the superscripts + and - refer to Ions and electrons respectively. 

Here n is the species number density, the average velocity. T the temperature, 

i 


I 

5 * 



I 




e the absolute value of electron charge, cp the electric field potential, m 
the mass, and t is the average collision time with the background gas. 

The upper sign choice in these equations is for ions, the lower is for electrons. 
In writing these equations we have assumed that the plasma is quasineutral 
(n + s» n = n) and the temperature of each species is maintained spatially 
uniform by frequent collisions with the background neutral gas. We also 
assumed that the frequencies of the motion we study are much less than the 
collision frequencies so that the inertia term can be ignored. The electric 
field is presumed electrostatic in nature. 

We can solve for nv± from Eq. (2.2) and substitute the results into 
Eq. (2.1). We obtain two equations for the determination of n and cp. In 
cylindrical coordinates these equations have the form 
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[ 




where V x H = 0 has been assumed. Here for each species b = eT/m is the 
mobility, D = Tt/id the diffusion coefficient, and b^ and represent b and D 
each divided by 1 + (Or ) 2 where fl = eH/mc is tne cyclotron frequency. The 
magnetic field is assumed to be so strong that (frr ) 2 » 1 for both electrons 
and ions. 

We shall repeat here the definition of the equilibrium state as given 
in I. Let us consider an equilibrium in which the density and the potential 
&re functions of r only. Eq. (2.3) reduces to 



where N and V are the equilibrium density and potential respectively. Making 
use of the fact D x ~ H 2 ~ r' 2 as is b ± and eliminating the b ± terms between 
the two equations implied above, one finds that 


li( r 3 — ) 

P*dr' r dr ; 


= 0 


and hence 


N = A + B r " 2 

The corresponding equation for V is 

d V _ c onst. 
dr r* 3 !? 

Since N is positive in the region of integration and since V must vanish 
at both ends, we must have 


V = 0 


(2.4) 


Adjusting the constants on N to fit the boundary conditions, one has 
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(2.5) 


It can be shown that a necessary condition for the instability is that 
Nr 2 decrease in the outward direction. Hence, to assure that we have such 
a configuration, we require that 

s « 

or expanding 


6s (R+df~ 
(R+d)*-R* 


£s x 2_d 

s'' R 


( 2 . 6 ) 


Henceforth, we can use the approximation 





(2.7) 


Let us now consider a density perturbation ni (r, z, t) and a potential 
perturbation V x (r, z, t) superimposed on the above equilibrium state. No 0 
dependence of the perturbation quantities is assumed since we expect the most 
unstable mode to be of the interchange type with no variation along the 
magnetic field. Substituting 

n = N(r) + n x (r, z, t) 
cp = V(r) + V x (r, z, t) 

into Eq. (2.3) and making use of the properties of the equilibrium solution, 
we arrive at the following equations governing the perturbations. 
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1( oV)^(dz^^ - + 

- | ^[rlN^IblHl) - ^[ (N+ni)b ±^j 
+ ^[(N+i^b^ 1 ]} = 0 


Since D ~ H 2 ~ r 2 as is b , we may write 
x i 


* «**»48* 

■ b ^('?^ : ^ r5N l7 1 ^ +N JP ) + (nV)bl i f ^(r 2 N)|| 1 
+ <a-'T J; )bi(i 2 ^:U 2 n I || 1 )- ^ (nj^)J = 0 

By virtue of the boundary conditions, we expect dnj/dr to be of order n x /d 
and similarly dVx/dr to be of order Vi/d. In Eqs. (2.8) r-derivatives may 
therefore be brought through r-dependent coefficients and allowed to act 
directly on dn x /dr, dV x /dr, nx(SVx/dr) and n^Vx/dz). The neglected terms 
(involving r-derivat*ves of the coefficients) are small compared to those 
retained in the parameter d/R or 6s/s. Hence using Eq. (2.7) one obtains 

— 1 - I)i( /"»* tlh + JT e 2^1 - k“n( 

3* " 1( ^ ^ 1 1 T 1755 +bxMt 3? + I? ) 

- h *irz~ b i[^ (n xl? ' + ^ in ils 1 )) 


(2.8) 


0 
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The coefficient of each term is a very weak function of r in the annular 
region and may be assumed constant. Defining x = r - R and neglecting small 
quantities of order d/R, we may write 


dnj 

at 


+ 3?n. 
DlI - 1 1 + 
ax 


— + Jlfc 2^1 

3z a ; " e hl^z 


- , 2>V ^v, 

+ b A 8 ( T~x + “1 ) 
ox 3z 


- H T T* + b *(^ <n lH 1) + ^ (n lH >J 

+ H(^ (n lH 1) - = 0 


( 2 . 9 ) 


Eqs. (2.9) are the approximate equations governing the perturbation quantities 
nj. and Vi. 

Let us now consider infinitesimal perturbations, when terms involving 
products of n x and Vi become negligible. The resulting linearized equations 
and the boundary conditions suggest perturbations of the form sin(nnx/d) 
exp(iu)t + ikz) where n is any positive integer and k is the axial wave number. 
Substituting perturbations of the above form into the linearized version of 
Equ- (2.9) and solving for the real and imaginary parts of u), we find 4 

- t [4 (T_b ^- T+b ^ ) ♦ If < d i- » i>] <b 1 + bi r 1 



(k 2 + X*f (Dxbl-t-Dlbi) -( 2&8/8deR) (kc/H)* (T*+T~) 
( **)( b++ bi) 


(2.10) 


where < = nn/d. The imaginary part of uj is negative at large value of H 
and the instability can occur. The neutrally stable state occurs when 


( nVun'r') 


(k*+ it ) (n»rf R a 

kV 2 1SS 


(2.11) 
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The right-hand, side is a minimum for k = X = ir/d. Hence the critical field 
H c at which the instability sets in is given by 

(BVHflPn = 


which is much larger than unity, in accord with our initial assumptions. 
Solving for H c , we find 



2 t T 2 C 2 R s 
b + b™ d "Ts 


( 2 . 12 ) 


In Fig. 2, we plot on the H-k plane the neutral stability curves representing 

</ = 0 for various values of n. The plasma is linearly stable for H < H . 

c 

Modes with small wavenumber do not produce sufficient polarization electric 

field for their growth while modes with too large a wavenumber are damped 

out by diffusion. So the instability first sets in only with the critical 

wavenumber k = rr/d and n = 1 at H = H . 

c c 

*>. Machine Calculation 

In this paper, we are concerned with the development of finite amplitude 
motion for H > H c « We shall integrate the full Eqs. (2.9) in time on a 
computer, assuming small initial, perturbations, and look into the ultimate 
fate of these perturbations. 

The perturbation equations will first be put into dimensionless form. 

We take the characteristic length in the x and z directions to be d and Ad 
respectively, where A is an adjustable parameter. The characteristic time 
is taken to be the approximate period t* of the critical mode in the limit 
of b + « b and T + » T 


> 
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-t* _ 2d eH 3 
T+ c 5s 


Making use of the following dimensionless variables, 

T1 = x /d 

§ = z/(Ad) 
t = t/t* 
p = nx/fis 
i = Vx/ (^- + ) 

We obtain from Eqs. ( 2*9 )> 

L£ _ A , jfp + J_ ifP) _ , xj£ 

3T Al '»(p A' if Az A 3f 

+ ^•)P ) 

+ A 3^(P|^ > + 

B l 1 }^ + A l l^ ) " Bz aI^ 

+ B 3 ( |^ + 


B,<P$>+^<Pff>] = 


where 


ov+H s 

2bj -c 13 


. 11 - 2 . 
‘ R is 


2 H ,+ 
~ b J- 
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R - 2H_s_/ y + T~ V 
B 1 ~ c 6s bj - T + ' 


“ 1 t£ ( 1 + 


T' 

T 


b 3 = f 3s<»: + »d 
b 4 = f (bi + bi) 


Eq. (3*3) is obtained from the ion version of Eqs. (2.9) and Eq. (3*^+) 
is obtained by eliminating dni/dt between the two Eqs. (2.9)* 

Let the density perturbation p(T], §, t) and the potential perturbation 
Y(T], 5, t) be expressed by 

pm, 5,t) = £ £ P nyn (T) sin (nr)?) exp(innr^) 

n=i (3.6) 

H'v 1*5 »T) = £ y I (T) sin(nTrl')) exp(iraTT-£) 

n~i h^-oo 

Here the choice of the sine series in the x direction satisfies the boundary 
conditions, and we use a discrete wave spectrum with a fundamental wave 
number of rr/(Ad) to approximate the continuous wave spectrum in the z direction. 

We substitute the Fourier series in (3*6) into Eqs. (3*3) and (3«^) and 
obtain the following system of equations 


- V 2 (n Z + < 1« + izn f T*. 

JTAjfnU n m+ 2 + iTT^W^ 
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B^[n Z * ) 2 ]p m + iE,rr fP nm 

♦ ¥ j (» z ’(aI j ]L 

+ 2 " B 4 A v n>J = 0 

(n = 1, 2, 3 , m = 0, ± 1, ± 2, ) 

where U , V and W are terms bilinear in p . . and Y., , and represent 
nm nm nm lj ij 

interactions between various modes. 



Here the summation on j is over odd integers if l ± n is even or zero and is 
over even integers if l ± n is odd. 


i 



Here sgn(n) = n/ In I if n ^ 0 and sgn(0) = 0. 

Computations restrict us to a finite sub-system of these equations. 
We use the finite sum 


P< 1 .?.r) 

if ( >i ,?,r) 


-e6 

n =l m=-M 

-t t 

n-i ms- M 


p (r) sin(n-rr>|) exp(imTT^) 


sin(n7D]) exp(imTrjf) 


(3.8' 


(3*9' 
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to approximate the series in (5-6) and restrict n and m in Eqs. (3-7) and 

(3*8) to 1 ^ n £ N and 0 £ j m | £ M. The series in the expressions for U , 

V and W are simply terminated by choosing p . = Y. = 0 whenever 

nm nm J ij ij 

i > N or I J I >M . We also make use of the reality conditions, p* = p , v 
• 1 ’ r nm K n(-m; 

and Y* = Y , \ where the asterisk indicates the complex conjugate. Thus 

nm n' v -m; ° 

we can further restrict the index m to 0 £ m £ M and obtain 2 X N X (2M + l) 

real equations from Eqs. (3 -T ) and (3-8) for 2 X N x (2M + l) real unknowns, 

namely, (p , p!\ p *) and (Y , Y Y * ) where l^n^N, 1 £ m £ M 
no nm nm no nm nm 

and the superscripts R and I denote the real and imaginary parts respectively. 

In carrying out the time integration, we used Hamming's modified predictor- 

corrector method with a special Runge-Kutta procedure for starting values. 

It is a stable fourth-order integration procedure. The routine for evaluating 

the time derivative 3 p^/dr as ^°H° WS * First, the spectral density 

p for all n and m are given as initial values or obtained from the previous 

time step. We then solve the system of linear algebraic Eqs. (3*8) in the 

unknown potentials Y • Second, d p /&r is then evaluated using Eq. (3 - 7 ) - 

nm nm 

It is essential that the solutions of Eqs. (3 -8 ) be sufficiently accurate 
in each step that rounding errors can never accumulate to any great extent. 

We used a Gaussian reduction method with iterative refinements. 5 The diagonal 
dominant property of the coefficient matrix enables us to obtain solutions 
with 6 digit accuracy in the norm of the solution vector in a few iterations. 

The time step used in integration must also be smaller than the period of 
the fastest motion of the system. We chose the time step such that the 
estimate of the local truncation error given by the predict or- correct or 
method is always less than 10 6 . It is also necessary that N and M be chosen 
sufficiently large that enough modes be included in the calculations. We 
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have no mathematical criterion for the requirement on N and M. However, one 
can see from Fig. 2 that modes with n > 2 are heavily damped when the 
instability first sets in for modes of n = 1. We therefore chose N = 2 
and neglected any modes with n > 2. In addition, one can see that, for 
H > H c , only modes with an axial wavenumber k in the band k < k 2 have 
positive growth rates. We chose our modes centered around these modes. 

Choosing M = A 2 with A = 2, 3, 4 successively, we increased the number of 
modes included as well as the fundamental wavelength. A finite number of modes 
might be sufficient to describe the system if modes at both ends of the 
spectrum remain at low level during the time of interest. 

The following numerical values of system parameters are used in computation, 
b /b + = 10 2 , T /T + = 20, 6s/s = 10 3 and d/R = 10 5 . Notice that the 
inequality (2.6) is satisfied. 

4. Results in the Slightly Super-Critical Regime 

We first investigated the evolution of an initial white-noise-like 

disturbance when the magnetic field is slightly above the critical value H c » 

We assumed as initial values Ip | = 10 3 for all n and m with the initial 

,K nm ! 

phase determined by random numbers. 

Fig. 3 shows the time development of various modes at H = 1.05 H 

c 

under the assumption of N = 2, A = 2 and M = 4. The numbers of each curve 
represent the mode number n and m defined in Eq. (3*9)» Initially the various 
modes are observed to grow or decay exponentially in time. This is what one 
expects since the initial perturbations are small and linear analysis giving 
exponential growth or decay should be valid. The initial growth rate and 
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oscillating frequency of each mode are in good agreement with that given by 
Eq. (2.10) of the linear analysis. Thus we have a good check on our computing 
codes. One notices in Fig. 3 the growth of the mode with n = 2, m = 0 as 
well as the critical mode n = 1, m = 2. The mode (2, 0) represents a modification 
to the equilibrium density. Its sign (not shown in the figure) is negative 
so as to reduce the equilibrium density in the region near the inner cylinder 
and increase it near the outer wall. That is, it acts to flatten the density 
gradient which is what we expect the nonlinear correction to do. Due to this 
change in equilibrium density, the critical mode (1, 2) levels off subsequently 
and a final steadily oscillating state with the saturation amplitude shown 
is reached. The final oscillating frequency is also shifted from that given 
by linear analysis evaluated at the critical field. 

We ran another computation under the same conditions, however, with the 
number of modes increased by choosing N = 2, A = 3* M = 9* The results are 
shown in Fig. 4. The same final oscillating state as in the previous case 
is reached with the critical mode labeled (1,3) because of the different 
value of A assumed. 

Fig. 5 shows the results obtained at H = 1.1 H c * Notice that three 
modes (l, 2), (l, 3)» and (l, 4) are located in the unstable wave band and 
they all have initial positive growth rates. Some of the modes like (l, 5 ) > 

(1, l), and (2, l) which decay initially are excited through mode interactions. 
However all modes finally subside except the oscillating critical mode (l, 3) 
with an axial wavenumber k = n/d and the time independent modification given 
by the mode (2, 0). 
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We also ran computations at other field values. To summarize our 
computed results in the slightly super-critical regime, we plot as functions 
of magnetic field the difference 6uj between the final oscillation frequency 
and the linear oscillation frequency evaluated at the critical field (Fig. 6), 
the squared final oscillation amplitude of the density wave and the magnitude 
of the time independent change in the equilibrium density (Fig. 7), and the 
squared final oscillation amplitude of the potential wave and the magnitude 
of the time independent change in the equilibrium potential (Fig. 8). In 
these figures, we also plot the analytic results obtained from Eqs (3. 10), 
(3-11), (4.7)> (4.8), and (2.23) of I. The analytic results are valid in 
the limit A = (H-H^)/^ « 1. It is seen indeed that the analytic results 
and our numerical results are in good agreement in this limit. 

5 . Results at Higher Magnetic Fields 

We tried next to extend our computation with the magnetic field increased 
further. To include more modes in the unstable wave band kj. < k < k 2 , 
we chose A = 4 and M = 16. Fig. 9 shows the time development of various modes 
at H = 1.3 H • Only modes with amplitude larger than 10 4 at t = >0 are 
shown in the figure. It is seen that initially each mode grows or decays 
much more rapidly because of the increased magnetic field. Some of the modes 
which decay initially begin to grow in a short time. These include uhe mode 
(2, l) and some other modes not shown in the figure. After a period of 
interactions, all these modes as well as those which grow initially begin to 
subside leaving a steady oscillation of the mode (1, 3 ) with an axial 
wavenumber k = 1.25 and the mode (2, 0) giving the time independent 

change of the equilibrium distribution. Fig. 10 shows a similar evolution 


19 


of modes at H = 1.6 H . In this case, however, the final oscillation is 

c 

dominated by mode (l, 6) which has an axial wavenumber k = 1.5 n /d. The 
modes (1, 5 ) and (l, 6) happen to be the modes with the largest growth rates 
at H = 1.3 H c and H = 1.6 H c respectively. Their rapid growth seems to 
enable them to take over the other modes. 

To study more closely the problem of wavenumber selection, we ran some 

computations under the assumption of a weighted initial disturbance. That 

is, we put the initial magnitude of one of the modes in the unstable wave 

band kj. < k < k 2 at larger magnitude than the rest. To be more specific, we 

put the weighted amplitude at 2 x 10 3 and the other amplitudes at 10 4 . 

Fig. 11 shows results obtained at H = 1.3 H with the assumption of such a 

c 

weighted initial amplitude for the mode (l, 6). It is seen that stable final 
oscillation with the dominance of the mode (l, 6), instead of the mode (l, 5)» 
can be established. The weighted initial amplitude of the mode (1, 6) is 
sufficiently strong that it can dominate and force down other modes which have 
larger initial growth rates. Fig. 12 shows the results obtained at H = 1.3 H c 
with the mode (1, 7 ) weighted. Contrary to the previous case, final oscillation 
with the mode (l, 7 ) cannot be reached. The mode (l, 7 seems unstable with 
respect to perturbation of other modes and the mode (l, which has the 
largest linear growth rate, takes over and finally dominates. Out of the 
seven modes (m = 2, 3, — 8) inside the unstable wave band at H - 1.3 H^, 
we find it possible to establish final stable oscillations of the mode with 
m = 4, 5* or 6 if the particular mode has its initial amplitude weighted 
with respect to the rest as described. Similarly, we find that, at H = 1.6 H c , 
it is possible to establish stable final oscillation of the mode with 
m = 4, 6, 7, or 8 out of the ten modes (m = 2, 3* H) inside the 


unstable wave band. 
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6 ' Discussion 

We would like to point out that the final state of the plasma indicated 
by our computed results is in sharp contrast with that obtained by Sato and 
Tsuda 2 for the case of the cross-field instability. Their computed results 
indicate that the instability develops explosively into a strong turbulence, 
while we hsve the dominance of a single finite amplitude mode up to H = 1.6 H c 
We do not have direct experimental sup- ort for the dominance of the single 
mode. However recent experimental work 6 on a similar macroscopic instability 
in a weakly ionized plasma, namely the spiral instability of the positive 
column, indicates the dominance of a single mode for magnetic field strengths 
up tc many times the critical value. Transition from laminar convection to 
turbulent convection could possibly occur at higher magnetic fields than those 
we investigated here. We do not go to higher magnetic fields because of the 
increased number of modes needed and the increased computing time. 

We would also like to mention that the idea of studying the wavenumber 
selection problem using a computational approach has been mentioned in the 
work of DiPrima and Rogers. 7 The nonuniqueness of the wavelength of the 
supercritical flow also occurs in some nonlinear hydrodynamic instabilities. 
The experimental work of Snyder 8 indicates the possibility of obtaining 
Taylor-vortex flows of different wavelengths. 

In summary, our numerical investigation shows that a weakly ionized 
plasma in a strong and nonunifora magnetic field is subject to an instability. 
Upon the onset of the instability, convective flow develops in the plasma 
and a final steadily oscillating state, dominated by a single finite amplitude 
mode together with a time independent modification of the original equilibrium 
is reached. The final oscillating state approaches that given by the 
analytic results of Simon 4 in the limit H -• H . We find the dominance of a 
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single mode up to H = 1.6 H c - However, the wavelength of the final oscillation 
reached is nununique and depends on initial conditions. For initial white- 
noise-like disturbances, we find that the mode with the largest linear growth 
rate will force down other modes and dominate. For weighted initial 
disturbances, we find there exists a subinterval of wavenumbers inside the 
unstable wave band. Final stable oscillation with a certain wavelength in 
this subinterval can be established, if the initial disturbance has a 
sufficiently strong component at the particular wavelength. 
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Figure 1. 
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Figure k. 


Figure 5* 


Geometry of the problem and picture of instability mechanism. 
Darkened areas indicate density enhancement due to perturbation. 
+ and - signs indicate polarization charges resulting from 
oppositely directed drifts of electrons and ions. 

Neutral stability curve of magnetic field H versus axial 
wavenumber k for various radial wavenumbers n. The plasma 
is stable for H < H c « The dashed line represents a band of 
unstable wavenumber at the operating magnetic field. 

Evolution of Fourier components of density wave at E = 1.05 H c 
for the case of white-noise-like initial disturbance. The 
curves are labeled according to the mode numbers n and m. 

Those modes not shown in the figure decay to less than 10 4 
before t = 1. 

Evolution of Fourier components of density wave at H = 1.05 H c 
for the case of white-ncise-like initial disturbance. The 
curves are labeled according to the mode numbers n and m. 

Those modes not shown in the figure decay to less than 10 4 
before t = 2 . 

Evolution of Fourier components of density wave at H = 1.1 H c 
for the case of white-noise-like initial disturbance. The 
curves are labeled according to the mode numbers n and m. 

Those modes not shown in the figure decay to less than 10 4 


before t = 3 


Figure 6. 


Figure 7. 


Figure 8. 


Figure 9. 


Figure 10. 


Comparison of the frequency shift 5ui as a function of magnetic 
field H from analytic results (line) and that from present 
numerical results (circles). 600 is the difference between 
the final oscillation frequency and the linear oscillation 
frequency evaluated at the critical field u> c . 

Comparison of the squared final amplitude of the density 

wave |piyj 2 and the magnitude of the time independent modification 

to the equilibrium density |p 2 o| as functions of magnetic 

field H obtained from analytic results (lines) and from 

present numerical results (squares and circles). 

Comparison of the squared final amplitude of the potential 
wave lY^I 2 and the magnitude of the time independent 
modification to the equilibrium potential |Y 2o | as functions 
of magnetic field H obtained from analytic results (lines) 
and from present numerical results (squares and circles). 

Evolution of Fourier components of density wave at H = 1.3 H c 
for the case of white-noise-like disturbance. The curves are 
labeled according to the mode numbers n and m. Those modes 
not shown in the figure have magnitude less than 10 4 for 
t * 30. 

Evolution of Fourier components of density wave at H = 1.6 H c 
for the case of white-noise-like disturbance . The curves 
are labeled according to the mode numbers n and m. Those 
modes not shown in the figure have magnitude less than 10 4 


for t i 40. 


Figure 11. 


Evolution of Fourier components of density wave at H = 1.3 
for the case of weighted initial disturbance. The curves are 
labeled according to the mode numbers n and m. The mode 
with n = 1 and m = 6 is the weighted component. Only the mode 
(2, 0) and modes with initial positive growth rate are shown 
in the figure. 

Figure 12. Evolution of Fourier components of density wave at H = 1.3 H c 
for the case of weighted initial disturbance. The curves 
are labeled according to the mode numbers n and m. The mode 
with n = 1 and m = 7 is the weighted component. Only the 
mode (2, 0) and modes with initial positive growth rate are 


shown in the figure. 































